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(57) Abstract: A method is provided for scaling up penneabilities assodated with a fine;-scale grid of cells representative of a porous 
medium (fig. 1) to permeabilities associated with an unstructured coarse-scale grid of cells representative of the porous medium (fig.3). 
An arcally unstructured, Voronoi, computational grid is generated using the coarse-scale grid as the genesis of the co^^)utational 
grid. The computational grid is then populated with permeabilities associated with the fine-scale grid. Flow equations are developed 
for the computational grid, the flow equations are solved, and inter-node fluxes and pressure gradients are then computed for the 
computational grid. These inter-node fluxes and pressure gradients are used to calculate inter-node average fluxes and average pres- 
sure gradients assodated with the coarse-scale grid. The inter-node average fluxes and average {aessure gradients associated widi 
the coarse grid are then used to calculate upscaled peraaeabilities associated with the coarse-scale grid. 
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METHOD OF UPSCALING PERMEABILITY FOR 
UNSTRUCTURED GRIDS 



FIELD OF THE INVENTION 

This invention relates generally to simulating fluid flow in a porous medium 
5 and, more specifically, to a flow-based method of scaling-up permeability associated 
with a fine-grid system representative of the porous medium to permeability 
associated with a coarse-grid system also representative of the porous medium. 

BACKGROUND OF THE INVENTION 

Numerical simulation is widely used in industrial fields as a method of 
10 simulating a physical system by using a computer. In most cases, there is a desire to 
model the transport processes occurring in the physical systems. What is being 
transported is typically mass, energy, momentum, or some combination thereof By 
using numerical simulation, it is possible to reproduce and observe a physical 
phenomenon and to determine design parameters without actual laboratory 
1 S experiments or field tests. It can be expected therefore that design time and cost can 
be reduced considerably. 

One type of simulation of great interest is a process of inferring the behavior 
of a real hydrocarbon-bearing reservoir firom the performance of a numerical model 
of that reservoir. The objective of reservoir simulation is to understand the complex 

20 chemical, physical, and fluid flow processes occurring in the reservoir sufficiently 
well to predict future behavior of the reservoir to maximize hydrocarbon recovery. 
Reservoir simulation often refers to the hydrodynamics of flow within a reservoir, but 
in a larger sense reservoir simulation can also refer to the total petroleum system 
which includes the reservoir, injection wells, production wells, surface flowlines, and 

25 surface processing facilities. 

The principle of numerical simulation is to numerically solve equations 
describing a physical phenomenon by a computer. Such equations are generally 



wo 00/79423 



PCTAJSOO/17101 



•2- 

ordinary differential equations and partial differential equations. These equations are 
typically solved using numerical methods such as the finite difference method, the 
finite element method, and the finite volume method among others. In each of these 
methods, the physical system to be modeled is divided into smaller cells (a set of 
5 which is called a grid or mesh), and the state variables continuously changing in each 
cell are represented by sets of values for each cell. An original differential equation is 
replaced by a set of algebraic equations to express the fundamental principles of 
conservation of mass, energy, and/or momentum within each smaller unit or cells and 
of mass, energy, and/or momentum transfer between cells. These equations can 

1 0 number in the millions. Such replacement of continuously changing values by a 
finite number of values for each cell is called "discretization". In order to analyze a 
phenomenon changing in time, it is necessary to calculate physical quantities at 
discrete intervals of time called timesteps, irrespective of the continuously changing 
conditions as a function of time. Time-dependent modeling of the transport processes 

1 S proceeds in a sequence of timesteps. 

In a typical simulation of a reservoir, solution of the primary unknowns, 
typically pressure and phase saturation or composition, are sought at specific points in 
the domain of interest. Such points are called "gridnodes" or more commonly 
"nodes." Cells are constructed around such nodes, and a grid is defined as a group of 

20 such cells. The properties such as porosity and permeability are assumed to be 
constant inside a cell. Other variables such as pressure and phase saturation are 
specified at the nodes, A link between two nodes is called a "connection." Fluid flow 
between two nodes is typically modeled as flow along the connection between them. 
In conventional reservoir simulation, most grid systems are structured. That 

25 is, the cells have similar shape and the same niunber of sides or faces. Most 

conmionly used structured grids are Cartesian or radial in which each cell has four 
sides in two dimensions or six faces in three dimensions. While structured grids are 
easy to use, they lack flexibility in adapting to changes in reservoir and well geometry 
and often can not effectively handle the spatial variation of physical properties of 

30 rock and fluids in the reservoir. Flexible grids have been proposed for use m such 
situations where structured grids are not as effective. A grid is called flexible or 
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unstructured when it is made up of polygons (polyhedra in three dimensions) having 
shapes, sizes, and number of sides or faces that can vary from place to place. 
Unstructured grids can conform to complex reservoir features more easily than 
structured grids and for this reason xmstructured grids have been proposed for use in 
5 reservoir modeling. 

One type of flexible grid that can be used in reservoir modeling is the Voronoi 
grid. A Voronoi cell is defined as the region of space that is closer to its node than to 
any other node, and a Voronoi grid is made of such cells. Each cell is associated with 
a node and a series of neighboring cells. The Voronoi grid is locally orthogonal in a 

10 geometrical sense; that is, the cell boundaries are normal to lines joining the nodes on 
the two sides of each boundary. For this reason, Voronoi grids are also called 
perpendicular bisection (PEBI) grids. A rectangular grid block (Cartesian grid) is a 
special case of the Voronoi grid. The PEBI grid has the flexibility to represent 
widely varying reservoir geometry, because the location of nodes can be chosen 

IS freely. PEBI grids are generated by assigning node locations in a given domain and 
then generating cell boundaries in a way such that each cell contains all the points 
that are closer to its node location than to any other node location. Since the inter- 
node cotmections in a PEBI grid are perpendicularly bisected by the cell boundaries, 
this simplifies the solution of flow equations significantly. For a more detailed 

20 description of PEBI grid generation, see Palagi, C.L. and Aziz, K.: "Use of Voronoi 
Grid in Reservoir Simulation," paper SPE 22889 presented at the 66th Aimual 
Technical Conference and Exhibition, Dallas, TX (Oct. 6-9, 1991). 

The mesh formed by connecting adjacent nodes of PEBI cells is commonly 
called a Delaunay mesh if formed by triangles only. In a two-dimensional Delaxmay 

25 mesh, the reservoir is divided into triangles with the gridnodes at the vertices of the 
triangles such that the triangles fill the reservoir. Such triangulation is Delaunay 
when a circle passing through the vertices of a triangle (the circiuncenter) does not 
contain any other node inside it. In three-dimensions, the reservoir region is 
decomposed into tetrahedra such that the reservoir volume is completely filled. Such 

30 a triangulation is a Delaunay mesh when a sphere passing through the vertices of the 
tetrahedron (the circimisphere) does not contain any other node. Delaimay 
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triangulation techniques are well known; see for example U.S. patent 5,886,702 to 
Migdal et aL 

Through advanced reservoir characterization techniques, it is common to 
model the geologic structure and stratigraphy of a reservoir with millions of grid 
5 cells, each populated wdth a reservoir property that includes, but is not limited to, rock 
type, porosity, permeability, initial interstitial fluid saturation, and relative 
permeability and capillary pressure functions. However, reservoir simulations are 
typically performed with far fewer grid cells. The direct use of fine-grid models for 
reservoir simulation is not generally feasible because their fine level of detail places 

10 prohibitive demands on computational resources. Therefore, a method is needed to 
transform or to scale up the fine-grid geologic reservoir model to a coarse-grid 
simulation model while preserving, as much as possible, the fluid flow characteristics 
of the fine-grid model. 

One key fluid flow property for reservoir simulation is permeability. 

1 S Permeability is the ability of a rock to transmit fluids through interconnected pores in 
the rock. It can vary substantially within a hydrocarbon-bearing reservoir. Typically, 
permeabilities are generated for fine-scale models (geologic models) using data from 
well core samples. For simulation cells, the tieterogeneities of the geologic model are 
accounted for by determining an effective permeability. An effective permeability of 

20 a heterogeneous medium is typically defined as the permeability of an equivalent 
homogeneous medium that, for the same boundary conditions, would give the same 
flux (amount of fluid flow across a given area per unit time). Determining an 
effective permeability, commonly called permeability upscaling, is not 
straightforward. The main difficulty lies in the interdependent influences of 

25 permeability heterogeneities in the reservoir and the applied boundary conditions. 

Many diflFerent upscaling techniques have been proposed. Most of these 
techniques can be characterized as (1) direct methods or (2) flow-based methods. 
Examples of direct methods are simple averaging of various kinds (e.g., arithmetic, 
geometric and harmonic averaging) and successive renormalization. The flow-based 

30 techniques involve the solution of flow equations and account for spatial distribution * 
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of permeability. In general, the flow-based methods require more computational 
effort but are more accurate than the direct methods. 

An overview of different upscaling techniques is provided in the following 
papers: Wen, X. H. and Gomez-Hernandez, J. J,, "Upscaling Hydraulic 
5 Conductivities in Heterogeneous Media: An Overview," Journal of Hydrology^ Vol. 
183 (1996) 9t32; Begg, S.H.; Carter, R.R. and Dranfield, P., "Assigning Effective 
Values to Simulator Gridblock Parameters for Heterogeneous Reservoirs," SPE 
Reservoir Engineering (Nov. 1989) 455-465; Durlofsky, L. J., Behrens, R. A., and 
Bemath, A., "Scale Up of Three Dimensional Reservoir Descriptions," Paper SPE 

10 30709 presented at the Annual Technical Conference and Exhibition, Dallas, TX 
(Oct. 22-25, 1995); and Li, D„ Cullick, A., Lake, L. W., "Global Scale-up of 
Reservoir Model Permeability with Local Grid Refmement", Journal of Petroleum 
Science and Engineering, VoL 14 (1995) 1-13. The upscaling techniques proposed in 
the past were primarily focused on structured grids. A need exists for a method of 

1 5 upscaling permeabilities associated with a fine-scale geologic model to permeabilities 
associated with an imstructured, coarse-scale reservoir simulation model. 



SUMMARY 

A method is provided for scaling up permeabilities associated with a fine-scale 
grid of cells representative of a porous medium to permeabilities associated with an 

20 imstructured coarse-scale grid of cells representative of the porous medium. The first 
step is to generate an areally imstructured, Voronoi, computational grid using the 
coarse-scale grid as the genesis of the computational grid. The cells of the 
computational grid are smaller than the cells of the coarse-scale grid and each cell of 
the computational grid and the coarse-scale grid has a node. The computational grid 

25 is then populated with permeabilities associated with the fine-scale grid. Flow 

equations, preferably single-phase, steady-state pressure equations, are developed for 
the computational grid, the flow equations are solved, and inter-node fluxes and 
pressure gradients are then computed for the computational grid. These inter-node 
fluxes and pressure gradients are used to calculate inter-node average fluxes and 

r 
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average pressme gradients associated with the coarse-scale grid. The inter-node 
average fluxes and average pressure gradients associated with the coarse grid are then 
used to calculate upscaled permeabilities associated with the coarse-scale grid. 

In a preferred embodiment, the computational grid is constructed from the 
5 coarse-scale grid to produce inter-node connections of the computational grid that are 
parallel to the inter-node connections of the coarse-scale grid. The cells of the 
computational grid are preferably approximately the same size as the fine-scale ceils. 
The computational jgrid is preferably populated with permeabilities by assigning to a 
given node of the computational grid, a predetermined permeability of a cell of the 

10 fine-scale grid that would contain the given node's location if the computational grid 
were superimposed on the fine-scale grid. The flow equations that are developed for 
the computational grid are preferably single-phase, steady-state equations. The inter- 
node average fluxes and average pressure gradients associated with the coarse-scale 
grid are preferably calculated using only the inter-node connections of the 

15 computational grid that fall within a predetermined sub-domain of the computational 
grid, and more preferably such calculations are made using only the inter-node 
connections of the computational grid that are parallel to the inter-node coimection of 
the coarse-scale grid associated with the sub-domain. The permeabilities associated 
with inter-node connections of the coarse-scale grid are preferably determined by 

20 calculating the ratio of the inter-node average fluxes to the inter-node average 
pressure gradients that were calculated for the coarse-scale grid. 



BRIEF DESCRIPTION OF THE DRAWINGS 

The present invention and its advantages will be better understood by referring 
to the following detailed description and the following drawings in which like 
25 elements are given like niunerals and letters. 

Fig. 1 illustrates a two-dimensional, structured, fine-scale grid commonly 
used in geologic modeling of a porous medium. 

Fig. 2 illustrates a two-dimensional, unstructured, coarse-scale grid having 
36 cells of which four cells are shown having nodes designated I, J, K and L. 
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Fig. 3 illustrates a two-dimensional, unstructured computational grid 
generated from the coarse-scale grid of Fig. 2; the cells of the computational grid 
were generated using a grid-refinement ratio of 2. 

Fig. 4 illustrates 10 unstructured cells generated in the practice of this 
5 invention from three unstructured cells of Fig. 2 having nodes I, J, and K; such 10 
cells were generated from the three cells using a grid-refinement ratio of 3. 

Fig. 5 illustrates IS unstructured cells generated in the practice of this 
invention from the three imstructured cells of Fig. 2 having nodes I, J, and K; such 
IS cells were generated from the three cells using a grid-refmement ratio of 4. 
10 Fig. 6 illustrates a two-dimensional computational grid generated from the 

grid of Fig. 2 using a grid-refinement ratio of 4, and shows a diamond shaped sub- 
domain KILJ. 

Fig. 7 illustrates four coarse-grid cells of Fig. 2 having nodes I, J, K, and L 

and the diamond shaped sub-domain R of Fig. 6. 
IS Fig. 8 illustrates a computational grid that is a sub-domain of the grid in Fig. 

6 and only includes the region of integration R of Fig. 6. 

Fig. 9 illustrates a computational grid that is a sub-domain of the grid in 

Fig. 6, which includes a region of integration R of Fig. 6 and a buffer zone (skin 

region) around region R. 
20 The drawings illustrate specific embodiments of practicing the method of this 

invention. The drawings are not intended to exclude from the scope of the invention 

other embodiments that are the result of normal and expected modifications of the 

specific embodiments. 

DETAILED DESCRIPTION OF THE INVENTION 

2S This invention provides a flow-based method for upscaling permeabilities 

associated with a fine-scale, geologic grid system representative of a porous medium 
to estimated permeabilities associated with a coarse-scale, imstructured grid system 
representative of the same porous medium. The invention is particularly useful in 
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upscaling peimeabilities associated with a geologic model of a reservoir, typically 
represented by structured cells, to permeabilities associated with a coarse-scale, 
xmstructured, Voronoi grid (also known as perpendicxilar-bisection or PEBI grid). 
As a first step in carrying out the scale-up method, an areally, unstructured, 
5 Voronoi grid (often referred to herein as a "computational grid") is generated using 
the coarse grid as its genesis. In accordance with the practice of this invention, the 
inter-node connections of both the coarse grid and computational grid are aligned 
(parallel) with each other. The computational grid is then populated with 
permeabilities fi-om the fine-scale, geologic grid. The next step is to set up flow 

10 equations, preferably single-phase, steady-state pressure equations, for the 

unstructured computational grid, solve the flow equations, and compute inter-node 
fluxes and pressure gradients for the computational grid. These fluxes and pressure 
gradients are then used to calculate inter-node average fluxes and average pressure 
gradients associated with the coarse grid. Permeabilities associated with the coarse- 

IS grid connections are then calculated using the average fluxes and average pressure 
gradients calculated earlier for the coarse grid. 

One embodiment of the invention v^ll now be described with reference to the 
drawings. Although the drawings illustrate only two-dimensional (2-D) grid systems, 
this invention is not limited to 2-D grids. As discussed later herein, the invention can 

20 also be applied to three-dimensional (3-D) grids. 

Fig. 1 illustrates an example of a fme-scale, geologic grid for use in 
representing a porous medium such as an aquifer or hydrocarbon-bearing reservoir. 
While Fig. 1 shows an 8x8 Cartesian grid of 64 cells, it should be understood that a 
typical geologic grid could contain millions of cells. All of the cells in Fig. 1 are 

25 structured, which is typical of geologic grids, with each cell having a node at its 
center. In Fig. 1, and in the other drawings, the thick dots denote cell nodes; the 
continuous lines denote the boimdaries of cells; and in Figs. 1-5 and 7 the dashed lines 
between nodes are referred to as inter-node connections. The triangles formed by the 
dashed lines in these figures form Delaunay meshes. 

30 Fig. 2 illustrates an example of an unstructured, PEBI grid system containing 

coarse-scale cells suitable for use in reservoir simulation. The grids of Figs. 1 and 2 
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represent the same domain of a porous medium. In most reservoir simulation 
applications, the average size of the coarse-scale cells of Fig. 2 would be several times 
larger than the average size of fine-scale cells of Fig. 1 . A reservoir simulation model 
may for example involve from about 10 to 1000 geologic cells for each reservoir 
5 simulation cell, with the length scale changing from 10 meters or more in the geologic 
model to 100 meters or more in the reservoir sirnulation model. Although the cells 
depicted in Fig. 1 are smaller than the cells depicted in Fig, 2 (for sake of clarity in 
presentation) they could be even smaller in actual cases. Techniques for generating 
both fine-scale, geologic grids and coarse-scale simulation grids are well known. The 

10 inventors have discovered an effective flow-based upscaling method for transforming 
the permeabilities associated with the fine-scale grid (such as a structured grid of 
Fig. 1) to an unstructured, coarse-grid system (such as an unstructured grid of Fig. 2). 

Fig. 3 illustrates one example of a computational grid generated in the practice 
of this invention that can be used in upscaling the permeabilities fix)m a geologic 

IS model of Fig. 1 to an unstructured reservoir simulation model of Fig. 2. The first step 
in the method of this invention is to generate a computational grid for the same 
domain covered by Figs. 1 and 2. 

The computational grid is generated in such a way as to produce xmstructured, 
PEBI cells that are smaller than the coarse-grid cells. Each inter-node connection of 

20 the computational grid is parallel to an inter-node connection of the coarse grid. As 
will be discussed in more detail below with respect to the preferred embodiment, this 
parallel aligiunent facilitates the upscaling of pressure gradients and fluxes which is 
an important step in the method of this invention. v 

The computational grid is generated by assigning nodes in addition to nodes of 

25 the coarse-grid and regenerating PEBI cells. The additional nodes are assigned 

proportionally on all coarse-grid connections and within each Delaunay triangle. The 
same number of new nodes is added to each connection and any number of new nodes 
may be added. Increasing the number of nodes increases the number of smaller PEBI 
cells that can be formed from the coarse-grid cells. Increasing the number of smaller 

30 cells of the computational grid is referred to herein as increasing the refinement. The 
degree of refinement can be measured by counting the number of additional nodes per 
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coarse-grid connection. If one node is added per coarse-grid connection, the 
refinement ratio is 2 since it divides the coarse-grid connection into two equal 
minigrid connections; if 2 nodes are added per coarse-grid connection, the refinement 
ratio is 3, and if 3 nodes are added, the refinement ratio is 4, and so on. 
5 The PEBI grid of Fig. 3 was generated with a refinement ratio of 2. Referring 

to Fig. 3, nodes I, J, and K represent nodes of the three cells illustrated in Fig. 2 
having nodes I, J, and K. To generate the grid of Fig. 3, one additional node was 
introduced midway between nodes K and I, one additional node was introduced 
midway between nodes K and J, and one additional node was introduced midway 

10 between nodes I and J, and so on for all other inter-node connections of Fig. 2. The 
PEBI grid of Fig. 3 was then generated based on the original nodes of Fig. 2 and the 
newly added nodes. 

Examples of computational grids having refinement ratios of 3 and 4 are 
shown in Figs. 4 and S, respectively. Fig. 4 shows 10 PEBI cells generated fi:om the 

1 S three cells of Fig. 2 having nodes I, J, and K, representing a refmement ratio of 3 . To 
generate the PEBI cells of Fig. 4, two nodes 10 and 1 1 were added proportionately 
between nodes K and I; two nodes 13 and 14 were added proportionately between 
nodes K and J; two nodes 15 and 16 were added proportionately between I and J; and 
one node 12 was introduced vvithin the triangle formed by node K, I, and J. The 

20 Delaunay mesh (the dashed lines) formed by inter-node connections of the 10 cells of 
Fig. 4 consists of 9 equal-sized Delaunay triangles having equal angles. All 9 smaller 
Delaunay triangles of Fig. 4 are similar to the Delaunay triangle IJK of Fig. 2. 

Fig. 5 shows 15 PEBI cells generated fi-om the three cells of Fig. 2 having' 
nodes I, J, and K, representing a refmement ratio of 4. To generate the 15 PEBI cells 

25 shown in Fig. 5, three nodes 20, 21, and 22 were added proportionately between 
nodes K and I of Fig. 2; three nodes 23, 24, and 25 were added proportionately 
between nodes K and J; three nodes 26, 27, and 28 were added proportionately 
between I and J; and three nodes 29, 30, and 31 were proportionately added within the 
triangle formed by node K, I, and J. Nodes 29, 30, and 31 lie at the intersections of 

30 lines drawn parallel to sides of triangle KIJ through nodes 20-28. Hence, all inter- 
node connections of the PEBI cells in Fig. 5, are parallel to one of the inter-node 
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connections between nodes K and I, between nodes K and J, and nodes I and J of 
Fig. 2. For example, referring to Fig. 5, the inter-node connection between nodes 20 
and 23 is parallel to the inter-node connection between I and J of Fig. 2; the inter- 
node connection between nodes 23 and 24 is parallel to the inter-node connection 
5 between nodes K and J of Fig. 2; and the inter-node connection between nodes 29 and 
30 is parallel to the inter-node connection between nodes K and I of Fig. 2. 

For a given refinement ratio, the number of similar, smaller, Delaunay 
triangles of a computational grid that can be produced is n^ times the number of 
Delaunay triangles of the coarse-scale grid, where n is a desired integer refinement 

10 ratio. For example, a computational grid with a refinement ratio of 4 will have 

sixteen (4^) Delaimay triangles (the Delaunay grid of Fig. 5) that are similar to, and 
smaller than, the Delaunay triangle KIJ of Fig. 2. 

The refinement ratio can be any integer greater than 1 . For most reservoir 
simulation applications, as a non-limiting example, the refinement ratio may range . 

1 5 fi-om 2 to about 1 0, and more typically the ratio will range from 2 to 5. The desired 
refinement ratio will depend on the relative sizes of the geologic grid and the coarse 
grid and the desired number of computational grid cells. The refinement ratio is 
preferably chosen to generate computational grid cells of approximately the same size 
as cells of the geologic grid for the same domain of the reservoir, and more preferably 

20 the refinement ratio chosen produces computational grid cells slightly smaller than the 
geologic cells. 

The method of this invention generates PEBI cells when all angles of each 

o 

Delaunay-mesh triangle form coarse-grid connections that are smaller than 90 (acute 

0 

triangle) or when one angle is at most equal to 90 (right triangle). However, if any 

o 

25 angle of a Delaunay triangle is greater than 90 (obtuse triangle), a true perpendicular 
bisector does not exist. For such triangles, the grid-generation technique of this 
invention v^U produce connections that are not parallel to the overlying coarse-grid 
connections. If any such non-parallel connections are generated, they are preferably 
not used in upscaling permeabilities firom the computational grid to the coarse-grid as 

30 discussed in more detail below. The coarse-grid is preferably generated such that 
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there are no obtuse triangles in the Delaunay-mesh. This can be achieved by 
adjusting node locations in most situations. 

The minigrid-generation technique of this invention is not limited to a 
triangular connection mesh (Delaunay-mesh); it is also applicable to structured grids 
5 in which the connection mesh comprises rectangles. The present computational grid 
generation technique is also applicable to a connection mesh that comprises a 
combination of triangles and rectangles. 
Permeability Population 

Once the computational grid has been generated, the next step in the method 

10 of this invention is to populate penneabilities onto the computational grid. This 

populating step is carried out using predetermined permeabilities associated with the 
fine-scale, geologic model (Fig. 1 for example). Various populating approaches can 
be used. One approach assigns permeabilities fix)m the geologic grid to nodes of the 
computational grid as follows. For a given node on the computational grid, the 

1 S geologic grid is searched until the cell of the geologic grid is located that contains the 
node of the computational grid, assuming that the computational grid is superimposed 
on the geologic grid. The predetermined permeability of the geologic cell that • 
contains a given node of the computational grid is assigned as the permeability of that 
given node. Each inter-node connection permeability of the computational grid is 

20 then computed by averaging (preferably by harmonic averaging) the two node 

permeabilities forming the inter-node connection. Another approach to populating the 
computational grid with permeabilities is to assign to the midpoints of the 
computational grid's inter-node connections, permeabilities obtained from the 
geologic cells. For a given midpoint of the computational grid's inter-node 

25 connection, a cell of the geologic grid is located that contains the given mid-point 
(assuming that the computational grid is superimposed on the geologic grid) and the 
permeability of that geologic cell is assigned to the given midpoint. If the geologic 
model has a diagonal permeability tensor, the two areal permeabilities can be 
combined based on the connection direction, before population. 
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Fluid Flow Equations 

After the computational grid has been populated with permeabilities, the next 
step in the practice of this invention is to develop flow equations for the 
computational grid, preferably single-phase, steady-state pressure equations for each 
5 cell of the computational grid, and using an assumed set of boundary conditions, solve 
the equations for each computational grid cell. The single-phase, steady-state 
pressure equation is given as: 

V.|VP = 0 (1) 

10 where P is pressure and ^ is the permeability tensor. 

The following description makes use of mathematical symbols, many of which 

are defined as they occur throughout the text. Additionally, for purposes of 

completeness, a symbols table containing definitions of symbols used herein is 

presented following the detailed description. 
IS Assuming, but not limited to, a scalar connection permeability, the discretized 

form of Eq. (1) for PEBI grids is: 

20 where T is transmissibility, subscript j refers to the node of interest and subscript £ refers 
to all of its neighbors. The term transmissibility as used in this description refers to a 
measure of the capability of a given viscosity fluid to move across a cell boundary (or 
inter-node connection) under a pressure drop. More specifically, transmissibility is 
known to those skilled in the art as a measure of the ability of a fluid to flow between two 

kA 

25 neighboring cells within a porous medium. Transmissibility is expressed as — , where k 

As 

is the effective permeability of the porous medium, A is the area of the boundary between 
the neighboring cells, and As is the average or characteristic distance that the fluid must 
travel in moving between the two cells. 
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One acceptable set of boundary conditions exerts a constant pressure gradient 
in the flow direction and assumes no flow in the transverse direction (no-transverse- 
' flow boundary conditions). As shown in Fig. 6, fluid flow is assumed to be injected 
or produced only in the direction of the arrows (from left to right). Another set of 
5 boundary conditions suitable for the practice of this invention is linear boundary 
conditions. In linear boundary conditions, a constant pressure gradient is applied in 
the flow direction along all boundaries. This can result in flow across all boundaries 
which can help in upscaling permeabilities in isolated areas. The foregoing boundary 
conditions should not be construed as limitation of the invention; other suitable 

10 boundary conditions can also be used. 

The pressure equation is solved in each of the two principal directions for 2-D 
applications (and in three principal directions for 3-D applications). Fluxes and 
pressure gradients are then computed for all inter-node connections of the 
computational grid. Persons skilled in the art of reservoir simulation would be 

IS familiar with development of suitable single-phase pressure equations for each cell of 
the computational grid, solving the pressure equations, and computing fluxes and 
pressure gradients for all inter-node connections of the computational grid. See for 
example, Verma, S., "A Control Volume Scheme for Flexible Grids in Reservoir 
Simulation," paper SPE 37999 presented at the Reservoir Simulation Symposium, 

20 Dallas TX (June 8-1 1, 1997). 
Upscaling 

Once the fluxes and pressure gradients computed are determined for the 
computational grid, the average fluxes and average pressure gradients associated with 
connections of the coarse grid are computed. These pressure gradients and fluxes are 

25 averaged over predetermined integration sub-domains associated with each coarse-grid 
cell, preferably sub-domains associated with each inter-node connection of the coarse 
grid. The ratio of the upscaled flux to upscaled pressure gradient then gives the upscaled 
permeability. This upscaled permeability can then be used to compute transmissibility. 
A mathematical basis for upscaling permeabilities from a computational grid 

30 to the coarse grid will now be provided. This description makes use of several 

mathematical symbols, some of which are defined as they occur in this description. 
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Additionally, for piuposes of completeness, a symbols table containing definitions of 
symbols used herein is presented following the detailed description. 

Darcy's law for 1-D single-phase flow in direction s is given by: 



Vids (3) 

Integrating Eq. (3) with respect to volume for computing an average flux over 
the volume of sub-domain R gives: 



(4) 



Assuming constant viscosity, the upscaled permeability (k ) for sub-domain 
R in direction s at the coarse-scale grid can now be defined as: 

IS Equation (S) can be expressed as a ratio of upscaled flux to upscaled pressure gradient 
by the following relationship, which was described in a paper by Rubin, Y. and 
Gomez-Hernandez, J. J., "A Stochastic Approach to the Problem of Upscaling of 
Conductivity in Disordered Media: Theory and Unconditional Nimierical 
Simulations," Water Resources Research, Vol. 26, No. 4, pages 691-701, April, 1990: 



(6) 

Approximating the integrals in Eq. (6) on the fine scale and substituting 
cotmection volume Vy(y^=A^)y the following expression for upscaled permeability 
25 (k;) is obtained: 
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i 



i 



(7) 



The connection transmissibility can now be computed as: 




(8) 



10 



15 



20 



Equation (7) can be applied using pressure gradients and fluxes computed 
from the pressure equation solution on the PEBI computational grid. An important 
step is to determine a sub-domain or region of integration associated with each 
coarse-grid connection over which the computational grid permeabilities are to be 
upscaled. 

Referring to Figs. 6 and 7, for upscaling areal connection permeabilities on 
PEBI grids, a suitable sub-domain is a diamond-shaped sub-domain R consisting of 
two triangles (KIJ and IJL) of Fig. 2. Other suitable sub-domain shapes can also be 
used. Equation (7) can now be used to compute upscaled permeability for each sub- 
domain using part or all of the computational grides inter-node connections within the 
sub-domain using a predetermined selection criterion. In a simple but efficient 
approach, the upscaled permeability uses only the inter-node connections of the 
computational grid within sub-domain R that are parallel to the coarse-grid's inter- 
node coimection representative of sub-domain R. In another embodiment, all of the 
computational grid's inter-node connections in sub-domain R could be used in the 
upscaling. 

Referring to Fig. 7, the computational grid inter-node connections in sub- 
domain R are used to obtain an upscaled permeability for the inter-node connection 
between nodes I and J. Since the pressure equation is solved in two principal 
directions in the areal domain, two upscaled permeabilities are obtained. If no- 
transverse-flow boundary conditions are used, the final upscaled permeability is 
selected based on the areal direction with the highest average pressure gradient for 
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each coarse-@id area! connection. Other techniques could be applied to upscale the 
permeabilities using combinations of the two pressure equation solutions. For 
example, when linear boundary conditions are used, a linear combination of the two 
pressure gradients can be obtained as follows: 

5 

ds {dsj, yds)2 (9) 

Constants a and b are determined such that the resulting pressure gradient is 
aligned with the connection direction for which the permeability is to be upscaled. 
1 0 The upscaled permeability now can be expressed as: 



k = '^' . . — — — 



dV 



dP 

ds 



2. 



dV 

(10) 



The upscaling technique of this invention can also be applied to 3-D structured 
1 5 as well as 3-D layered PEBI grids (also referred to by some as 2Vi-D PEBI grids). 
Extension to 3-D unstructured grids is also possible. The layered PEBI grids are 
unstructured areally and structured (layered) vertically. One way to create such grids 
is to project 2-D areal PEBI grids on geologic sequence surfaces. Areal connection 
permeabilities can be upscaled layer-by-layer or over multiple layers. For upscaling 
20 vertical (inter-layer) connection permeabilities, the flow equations would be solved 
with the pressure gradient in the vertical direction. Construction of layered 3-D grids 
is described by (1) Heinemann, Z. E., et al., "Modeling Reservoir Geometry With 
Irregular Grids," SPE Reservoir Engineering, May, 1991 and (2) Verma, S., et al., "A 
Control Volimie Scheme for Flexible Grids in Reservoir Simulation," SPE 37999, 
25 SPE Reservoir Simulation Symposium, Dallas, TX, June, 1997. The upscaling 
method of this invention can also be extended to full 3-D PEBI grids. 
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A veiy useful aspect of this invention is that the scale-up method can be 
applied to upscale permeabilities for an arbitrary number of coarse-grid connections 
simultaneously. The fluid flow equations can be solved on a computational grid 
corresponding to a large coarse-grid domain (scale-up domain), thus permitting a 
5 large number of permeabilities to be upscaled at once. For example, the fluid flow 
equations can be solved for the entire computational grid of Fig. 6 to upscale 
permeabilities for all the coarse grid inter-node connections shown in Fig. 2. 
Alternatively, the fluid flow equations can be solved on a smaller computational grid 
corresponding to one or a few coarse-grid connections. An example of a 

10 computational grid for upscaling one coarse-grid connection permeability is shown in 
Fig. 8. Here, the computational-grid coincides with the diamond-shaped region of 
integration R used in upscaling permeability for coarse-grid connection IJ. 

The flow equations can also be solved on a larger computational grid than the 
scale-up domain thus allowing a buffer or skin around the scale-up domain whether it 

1 S consists of single or multiple connections. An example of this is shown in Fig. 9 
where the computational grid includes the diamond-shaped region of mtegration R 
(KILJ) for cormection IJ and a bu£fer zone or skin region surrounding region R. 
These options offer flexibility in selecting scale-up domain size as well as 
computational-grid size for upscaling, which could be important for scale-up accuracy 

20 in certain applications where actual field boundary conditions are different from the 
scale-up boundary conditions. 

Although the mathematical description of the upscaling method of this 
invention has been presented with a single-phase, no gravity, steady-state, scalar 
permeability formulation, the method by no means is limited to these assumptions. 

25 Persons skilled in the art can extend the scale-up method of this invention to include 
multiphase flow, unsteady state flow, a Ml permeability tensor, and gravity using 
the information disclosed in this patent. 

Although not shown in the drawings, in practicing the present invention, the 
properties of the domain being simulated can optionally be visualized on 2-D and 

30 3-D imstructured grids using suitable computerized visualization equipment. The 
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visualization of the properties can be helpful in analyzing application of the 
invention. In 2-D, a connection property (such as flux, pressure gradient or 
connection permeability) can be represented by displaying a diamond-shaped region 
associated with each connection in color from a predefined color scale. The apexes 
5 of this diamond-shaped region include the two node locations forming the 

connection and the end-points of the perpendicular bisector of the connection. In 
3-D, the Delaimay triangulation can be represented by a ball and stick network. 
Node properties (such as pressure, saturation, composition or porosity) can be 
displayed by coloring the balls using a color scale. Similarly, the connection 

10 properties can be displayed by coloring the sticks. The color scale can display the 
magnitude of a property. For displaying a vector property (for example, pressure 
gradient or flux), an arrowhead can be used to represent the flow direction. 

The principle of the invention and the best mode contemplated for applying 
that principle have been described. It will be apparent to those skilled in the art that 

15 various changes may be made to the embodiments described above without 
departing from the spirit and scope of this invention as defined in the following 
claims. It is, therefore, to be understood that this invention is not limited to the 
specific details shown and described. 



20 

Symbols 

a, b == Constants, dimensionless 

2 

A - Area, ft 

25 k - Permeability, md 

- = EffectivePermeability of Entire Grid Domain, md 

r = Pressure, psi 

T = Transmissibility, md-ft 

u = Darcy Flux, ft/day 

30 F = Volume, ft^ 

AP = Pressure Drop, psi 

As = Connection Length, ft 

// = Viscosity, cp 
Subscripts 
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c = Coarse Grid 

/ = Computational grid 

/ = Connection Index 

// = Node Index 

5 R - Sub-domain of Integration 
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What is claimed is: 

1 . A method of scaling up permeabilities associated with a fine-scale grid of cells 
5 representative of a porous medium to permeabilities associated with an 

unstructured coarse-scale grid of cells representative of the porous medium, 
comprising the steps of: 

(a) generating an areally unstructured, Voronoi, computational grid using 
the coarse-scale grid as the genesis of the computational grid, the cells of 

10 the computational grid being smaller than the cells of the coarse-scale 

grid and each cell of said computational grid and said coarse-scale grid 
having a node; 

(b) populating the computational grid with permeabilities associated with 
the fine-scale grid; 

1 5 (c) developing flow equations for the computational grid, solving the flow 

equations, and computing inter-node fluxes and pressure gradients for 
the computational grid; 

(d) using the fluxes and pressure gradients computed in step (c) to calculate 
inter-node average fluxes and average pressure gradients associated with 

20 the coarse-scale grid; and 

(e) calculating upscaled permeabilities associated with the coarse-scale grid 
using the average fluxes and average pressure gradients calculated in 
step (d). 

25 2. The method of claim 1 wherein inter-node connections of the computational 
grid cells are parallel to inter-node connections of the coarse-scale grid. 



3. 

30 



The method of claim 1 wherein cells of the computational grid are 
approximately the same size as the cells of the fine-scale grid. 
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4. The method of claim 1 wherein cells of the computational grid are smaller 
than the cells of the fine-scale grid. 

5. The method of claim 1 wherein the permeabilities populated in step (b) are 
5 assigned to nodes of the computational grid, the permeability assigned to a 

given node of the computational grid corresponding to a predetermined 
permeability of a cell of the fine-scale grid that would contain the given node's 
location if the computational grid was superimposed on the fine-scale grid. 

10 6. The method of claim 5 further comprises calculating the permeability for a 
given inter-node connection of the computational grid by harmonically 
averaging the permeabilities at the two nodes forming the given inter-node 
connection. 

.15 7. The method of claim 1 wherein the permeabilities populated in step (b) are 
assigned to midpoints of a given inter-node connection of the coniputational 
grid, the permeability assigned to the given inter-node connection 
corresponding to a predetermined permeability of a cell of the fine-scale grid 
that would contain the mid-point of such inter-node connection if the 

20 computational grid was superimposed on the fine-scale grid. 

8. The method of claim 1 wherein the calculation of the inter-node average 

fluxes and average pressure gradients in step (d) is determined by using only 
the fluxes and pressure gradients computed for the inter-node connections of 
25 the computational grid that fall within a predetermined sub-domain of the 

computational node. 



30 



9. 



The method of claim 1 wherein the calculation of the average fluxes and 
average pressure gradients in step (d) for a given inter-node connection of the 
coarse-scale grid is determined using only fluxes and pressure gradients 
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computed in step (c) for the inter-node connections of the computational grid 
that are parallel to the given inter-node connection. 

10. The method of claim 1 wherein the permeabilities calculated in step (e) are 
5 determined for a given node of the coarse-scale grid. 

1 1 . The method of claim 1 wherein the permeabilities calculated in step (e) are 
determined for a given inter-node connection of the coarse-scale grid by 
calculating the ratio of the average flux to the average pressure gradient 

1 0 computed in step (d) for the given inter-node connection. 

12. The method of claim 1 wherein cells of the fine-scale grid are structured. 

13. The method of claim 1 wherein the coarse-scale grid is a PEBI grid. 

15 

14. The method of claim 1 wherein both the coarse-scale grid and the 
computational grid are PEBI grids. 

15. The method of claim 1 wherein inter-node connections of the coarse grid forms 
20 Delaunay triangles and the computational grid generated in step (a) contains 

similar, smaller Delaunay triangles equal in number to times the nimiber of 
Delaunay triangles of the large-scale grid, where « is a predetemiined integer 
refinement ratio used to generate the computational grid. 

25 1 6. The method of claim 1 wherein all cells are three-dimensional. 

17. The method of claim 16 wherein the coarse-scale grid and the computational 
grid are both unstructured areally and structured vertically. 
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1 8. The method of claim 1 further comprises determining inter-node connection 
transmissibilities of the coarse-scale grid using permeabilities calculated in 
step (e). 

5 1 9. The method of claim 1 wherein the flow equations of step (c) are single-phase 
and steady-state. 

20. A method for estimating pemieability of each cell of a first grid having a 

multiplicity of cells of a subterranean geologic domain using a predetermined 
10 permeability for each cell of a second grid representative of the domain, said 

second grid containing a larger number of cells than the first grid, the method 
comprising: 

(a) constructing an unstructured, third grid representative of the domain 
comprising approximately the same or greater number of cells than the 

IS second grid, each cell of the first, second, and third grids having a node 

and each link between two nodes of adjacent cells being an inter-node 
connection, substantially all of the inter-node connections of the third 
grid being parallel to the inter-node connections of the first grid; 

(b) for each node of the third grid, assigning a permeability corresponding to 
20 the permeability of a cell of the second grid that contains the node 

location of the third grid; 

(c) developing a single-phase, steady-state pressure equation for each cell of 
the third grid system; 

(d) solving the pressure equations and computing fluxes and pressure 
25 gradients for all inter-node connections of the third grid; 

(e) computing an estimated permeability for a given connection of the first 
grid using inter-node connections of the third grid; and 

(f) repeating step (e) for all connections of the first grid. 



30 21. 



The method of claim 20 fiirther comprises computing the permeability in 
step (e) by the additional steps of determining average fluxes and average 
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pressure gradients over sub-domains associated with a given grid inter-node 
connection of the first grid and calculating a ratio of the average flux to the 
average pressure gradient, thereby obtaining the permeability for the given 
inter-node connection of the first grid. ^ 

5 

22. The method of claim 20 wherein the permeability computation of step (e) uses 
only inter-node connections of the third grid that are parallel to inter-node 
connections of the first grid. 

10 23. The method of claim 20 wherein all grid cells are three-dimensional. 



24. The method of claim 23 wherein the second grid and the third grid are each 
imstructured areally and structured vertically. 

15 25. A method for estimating permeabilities associated with cells of a large-scale 
grid representative of fluid flow in a porous mediiim using predetennined 
permeabilities associated with cells of a small-scale grid also representative of 
fluid flow in the porous medium, each cell of the large-scale grid having a 
node and each node being linked to adjacent nodes to form inter-node 

20 connections and such connections foiming Delaimay triangles, comprising the 

steps of: 

(a) constructing a computational grid by dividing each Delaunay triangle of 
the large-scale grid into a multiplicity of similar, smaller Delaunay 
triangles, the sides of such smaller Delaunay triangles being inter-node 

25 connections of the computational grid and the inter-node connections of 

the large-scale grid and the computational grid being aligned with each 
other; 

(b) assigning permeabilities to the computational grid corresponding to the 
predetermined permeabilities of the small-scale grid; 

30 (c) developing a single-phase, steady-state pressure equation for each cell of 

the computational grid, solving the pressure equations, and computing 
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fluxes and pressure gradients for all inter-node connections of the 
computational grid; 

(d) using the fluxes and pressure gradients computed in step (c) to calculate 
an average flux and an average pressure gradient for each inter-node 

5 connection of the large-scale grid; and 

(e) calculating a permeability associated with a given inter-node connection 
of the large-scale grid using the average flux and average pressure 
gradient calculated in step (d). 



10 26. The method of claim 25 wherein the number of nodes of the computational 

grid are approximately the same as the number of cells of the small-scale grid. 

27. The method of claim 25 wherein the number of nodes of the computational 
grid are smaller than the number of cells of the small-scale grid. 

15 

28. The method of claim 25 wherein the number of similar, smaller Delaunay 
triangles is nf' times the number of Delaunay triangles of the large-scale grid, 
where n is a predetermined integer refinement ratio. 
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